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INTRODUCTION 

Over the years, our research into turbulence at Edin- 
burgh has concentrated on the application of renormaliza- 
tion methods to the prediction of the energy spectrum of 
isotropic turbulence. General discussions of this work will 
be found elsewhere (McComb 1990, 1995), while accounts 
of specific progress have been given previously in this con- 
ference series (McComb & Shanmugasundaram 1983, Mc- 
Comb, Filipiak, Roberts & Watt, 1991). 

From a practical point of view, the most promising de- 
velopment in this area is undoubtedly Renormalization 
Group or RG. If we work in the Fourier representation, 
in principle, this involves the progressive averaging out of 
high-wavenumber modes in bands, with reseating at each 
step, until a fixed point is reached. The result is, in effect, 
a 'subgrid model' for large-eddy simulation. 

RG has enjoyed its successes in other areas of statistical 
physics. However, its application to turbulence faces sev- 
eral technical difficulties, which have to be circumvented by 
uncontrolled approximations. Indeed, in view of the deter- 
ministic nature of the Navier-Stokes equations, it is clear 
that the operation of averaging out the high-wavenumber 
modes while keeping the low-wavenumber modes constant, 
cannot be done rigorously and in itself can only be an ap- 
proximation. 

With points like this in mind, we have recently adopted 
direct numerical simulation as a tool for probing the basic 
feasibility of using RG techniques to reduce the number of 
degrees of freedom requiring to be numerically simulated. 
In this paper, we present some of the first results of this 
approach. We begin by discussing the RG approach in 
detail. 

RENORMALIZATION GROUP THEORY 

Basic Equations 

Working in Fourier-wavevector (k) space and restricting 
our attention to turbulent velocity fields which are homo- 
geneous, isotropic and stationary, we may write the pair- 
correlation of velocities as 

{u^{k,t)u^{k' ,t')) = Q{k,t~ t')D^^{k)5{k - k'), (1) 



where Q{k, t — t') is the spectral density and the projector 
Dap{k) — da0+kaki3k~^ ariscs duc to the incomprcssibility 
condition. Thus, the energy spectrum E{k) = 4'Kk^Q{k) 
with Q{k) — Q{k,0) and the maximum cut-off wave- 
number, fco, is defined via the dissipation integral 



dk 2uok^E{k) 



dk 2uok^E(k), 



(2) 



where e is the dissipation rate, vo is the kinematic viscosity, 
and fco is of the same order of magnitude as the Kolmogorov 
dissipation wave-number. 

Renormalization Group Theory 

Taking our goal to be the calculation of the energy spec- 
trum E{k), our intermediate objective is to find an analyt- 
ical method of reducing the number of degrees of freedom 
(or Fourier modes), in order to make the numerical solu- 
tion of the equations of motion a practical proposition. Let 
us consider how this might be done by using RG. 

First, we divide up the velocity field at = fci as 
Ua{k,t) — u~{k,t) for < fc < fci and Ua{k,t) = u^{k,t) 
for fci < k < ko, where ki — {1 — rj)ko and the bandwidth 
parameter rj satisfies the condition < 77 < 1. Work- 
ing with the standard form of the solenoidal Navier-Stokes 
equation in fc-space, we may write the evolution of the low- 
k velocity field for < fc < fci as 



;(k,t) 



= M^^^ik) J d'j ^UpQ,t)u-,ik~],t) 

+ 2u-;^{i,t)u+{k-i,t)+u+{i,t)u+{k~i,t) 



.(3) 



and the evolution of the high-fc velocity field for the first 
shell, fci < fe < fco, as 
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Figure 1: Convergence of the Kolmogorov spectral constant 
a to the fixed points for several values of the bandwidth pa- 
rameter r\. 



Figure 2: Dependence on the bandwidth parameter 77 of 
the calculated values of the Kolmogorov spectral constant a 
based on equation (^). 



(j, t)?i+(k - j, t) + ^t;^(j, t)u+(k - j, t) 



,(4) 



where the superscripts + and — on MQ,/3-y(k) have the 
same significance as for Ma(k,t), and the symmetrized 
inertia! transfer operator Ma,^-,(k) = (2i)~^\kfjDa'^{\s) + 

In principle, the RG approach involves two stages: (i) 
Eliminate the high-fc modes, u""", which appear in equation 
for < fc < fci, by solving for the mean effect of the 
high-fc field. This results in an increment to the viscosity, 
i.e. ^ v\ = vq -\-&vq. (ii) Rescale the basic variables, so 
that the Navier-Stokes equation for < fc < fci looks like 
the original Navier-Stokes equation for < fc < fco- 

Although this procedure is appealingly simple and has 
a clear physical interpretation, it has not proved easy to 
put into practice in the turbulence problem. A typical 
approach is to eliminate all the high-fc effects in equation 
(§)) by substituting the solution of equation (^), directly 
into the modes in the u~ equation. However, prob- 
lems are then encountered because of the mode coupling 
between u~ and u^. Even if one succeeds in carrying out 
the first part, the further problem of averaging out the 
high-fc modes arises immediately, because u~ and u""" are 
not statistically independent. This problem was avoided 
by Foster, Nelson and Stephen (1977; hereafter referred to 
as FNS) in their pioneering study of stirred fluid motion, as 
they restricted their attention to stirring forces which were 
multivariate normal and excluded the effects of the turbu- 
lence cascade. Ifowever, it has been shown that the use of 
a 'filtered' average by FNS to eliminate the u~ equation is 
really an uncontrolled approximation (Eyink, 1994). 

Iterative-Averaging RG with Results 

Here, we follow the method of iterative averaging, which 
is based upon the derivation of a recurrence relation and, 
in principle, eliminating finite blocks of modes (i.e. high-fc 
modes) while maintaining the form invariance of the dy- 
namical equation. Apart from the work of FNS, elimina- 
tion procedures can be performed by 'conditional' averag- 
ing, first introduced by McComb (1982). Further details 
about the conditional average have been given elsewhere 
(McComb, Robert and Watt, 1992). The basic ansatz of a 
conditional average is that a small uncertainty ("l"", say) 
at the cutoff wavenumber will generate chaotic behaviour 
for the high-fc modes. Although the introduction of $~ 
has been accepted, mainly due to the chaotic nature of the 
Navier-Stokes equations, it might be interesting to see how 
'rapidly' chaotic behaviour develops from the given small 



$~ by numerical simulation. This aspect is one of our cur- 
rent tasks and the results will be reported in due course. 

The current result of the iterative-averaging calculation 
for the Navier-Stokes equations after first eliminating the 
high-fc effects is 



(k,t) 



= M. 



(k) / d j tt^(j,t)u^(k-j,i). 



(5) 



where v\ = vq^ Suo{k) and 



^TV[M-,^(k)A/+,(k - j)D^„(j)] 



- I'ojk - j| 



(6) 



Here, we consider space dimension d — 3. This result can 
be extended to further shells, and we have shown elsewhere 
(McComb and Watt, 1992) that a fixed point is reached 
under numerical iteration of the recursion relations (see 
also Figure |l|). In Figure 0, we show a calculation of the 
Kolmogorov constant a =1.60 ± 0.01 independent of the 
bandwidth of modes being eliminated for bandwidths in 
the range 0.25 < r? < 0.45, in agreement with experiment. 

NUMERICAL SIMULATIONS 

Two programmes of numerical simulation are being car- 
ried out — one at the University of Edinburgh in the 
United Kingdom, the other at the Swiss Federal Insti- 
tute of Technology, Lausanne. A large number of runs 
have already been carried out at Lausanne, and this paper 
presents some of the results obtained so far. 

The simulations themselves are very similar, while the 
computer systems on which they are run differ greatly. At 
Edinburgh, work is carried out on a parallel machine, the 
Cray T3D, while in Lausanne a parallel- vector machine, 
the NEC SX-4, is used. 

The simulations discussed in this paper were carried out 
at a resolution of 256"^ , requiring approximately 14 seconds 
of SX-4 time per time-step on a single processor. 

The general method of such simulations has been well 
established. We follow the work of Orszag for the con- 
struction of initial velocity fields (1969) and in the use 
of a pseudospectral method (1971). The time integration 
scheme is a second-order Runge-Kutta method and partial 
dealiasing is achieved by way of a random-shifting method 
(see, for example, Rogallo, 1981). 
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Characteristics of the simulation 
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Initial Conditions 

The simulations are started with an initial energy spec- 
trum of the form 



E{k,Q) = 16(2/7r)2ii2fc-5it4gxp[_2(fc/fcp)2] 



(7) 



where fcp is the location of the spectrum's maximum and 
uq is the required initial r.m.s. velocity. 

Forcing 

Stationary turbulence is obtained by use of a determin- 
istic forcing term 

■'"^ ' ' \ otherwise, ^ ' 

where e is the mean dissipation rate, and 



Eiit) 



E{k,t)dk. 



(9) 



There is no preferred direction in this forcing and the tur- 
bulence rapidly reaches a statistically isotropic and steady 
state. 

Statistics 

While our simulations are entirely conventional, we do 
not rely solely on the usual practice (as justified by 
isotropy) of averaging over shells in wavenumber space 
in order to obtain statistical quantities, but also generate 
many realizations in order to increase our sample size. 

The main characteristics of the simulation are reported 
in Table ^ where At is the time step, T is the integration 
time, Vo IS the molecular viscosity, fc/ is defined in (^), fco 
is the ultraviolet cut-off, e is the mean dissipation rate, R\ 
is the Reynolds number based on the Taylor microscale, L 
is the integral scale, A is the Taylor microscale, te is the 
turnover time and S3 and S4 are respectively the skewness 
and flatness of the velocity derivative. 

The equations have been integrated for more than 60 
turnover times and about 200 box-realizations of each com- 
ponent of the velocity field have been stored in a database. 
Since these box-realizations are separated by ~ te/4: they 
can be considered statistically independent for the middle- 
range-scales and the small-scales. 

RESULTS 

We wish to assess the freedom to carry out conditional 
averages of the type required by RG. In principle we may 
do this by extracting, from an ensemble of realizations of 
the velocity field 



X = {Xi"^ (fe, t) I Q = 1, 2, 3; t G [0, T]; 

< |fc| < ko;n= 1, ...,iV}, 



(10) 



two disjoint subensembles y and Z chosen such that, for a 
prescribed C > Oj 



|y('"'(fe,t) - z^"'\k,t)\ 



for all < Ifcl < kc 



]y('")(fe,t)P) 

;m = l,...,M;te[0,T], (11) 




Figure 3: (a) Relative energy error for kt — 10, kc = 15, 
C = 0.5 and a — 1. (b) A selected set of realizations showing 
strong fluctuations for k > 15. 



for each realization y'""^ € V and Z^""^ £ Z. We may 
then define the relative energy of the error 

{{u{k,t)-w{k,t)f) 



r{\k\) 



2(u(fc,t)2) 



(12) 



where u{k,t) G y and w{k, t) G Z. fit is important to note 
that the averages in the definition ( [l2[ ) are, in this context, 
subensemble averages defined on v an d Z and not ensem- 
ble averages on A".) In equation (|l2| ) and subsequently, 
we assume that the fields are statistically stationary and 
isotropic, therefore r depends only on |fc|. Since the two 
fields are very close when < |fe| < kc, ''(Ifej) will be much 
smaller than 1 in this interval, indicating that the fields 
are almost completely correlated. If the error between the 
fields grows in such a way that they become decorrelated, 
we will have r(|fc|) — > 1 as |fc| > kc increases. 

In practice, our 200 box-realizations are not sufficient 
for the above analysis and we shall describe how we have 
extracted, using a partial sampling technique, enough real- 
izations to compute the relative energy of the error defined 
by (0). 

In order to this, we have performed the following partial 
Fourier transform of one component of the velocity field 



Uc,{x,y,k) = — 



Ua{x, y, z)e^''^ dk, 



(13) 



then we have selected, for each box-realization, a set of 
realizations, say Ua{xi,yi,k), where the spacing 5x = 
\xi+i — Xi\ — \yi+i — yi\ is chosen such that the realiza- 
tions are (approximately) independent for the range of k 
we consider (if we consider only the scales such that k > kb, 
then 5a:: — 2Ti/kh). The union of all these realizations ob- 
tained for each of the box-realizations will constitute our 
ensemble X. The subensemble y is formed by choosing an 
arbitrary subensemble of X. To select the subensemble Z, 
we impose the condition 



|y('")(fc) _^('")(fc)|2 
2(|F('")(fe,t)|2) 
for all ki, < k < kc \ m = 1, 



...,M. 



(14) 



Note that the time dependence does not appear in the 
equations since all the box-realizations used to form the 
ensemble X are taken in the statistically steady regime. 
Figure ^a) shows the relative energy error 



{(U{k) - W{k)f) 



(15) 



where u ^ y and w £ Z for fcj, = 10, fee = 15, C, = 0.5 
and a = 1. The number of realizations M is 2533. Though 
the number of realizations is not large enough to have a 
smooth converged solution, one can see that the relaxation 
to a chaotic regime is indeed very fast. Figure |(b) shows 
a selected set of realizations for which one can observe that 
the constraint imposed for 10 < fc < 15 does not prevent 
strong fluctuations for k > 15. The convergence of r{k) is 
difficult to improve, due to the restriction on the number 
of realizations available for a given constraint. 

Another natural way in which the small-scale proper- 
ties of a conditional subensemble may be investigated is 
by studying the probability density functions (pdfs) of ve- 
locity increments. In physical-space, we can use homo- 
geneity in the three dimensions and have sufficiently large 
subensembles to compute high-order statistics and pdfs. 
The velocity increments are defined by the following rela- 
tion 



Su{x, h) — u{x + h) — u{x), 



(16) 



where /i is a displacement vector and x the position. Since 
the fields are statistically isotropic, we can restrict our- 
selves to the study of the longitudinal velocity increment 
SvL{h) which is the projection of Su{h) on the direction 
of the vector h and the lateral velocity increment Syrih) 
which is the projection of 5u{h) on a direction perpendic- 
ular to h. For the purpose of this paper, we have only 
studied the longitudinal velocity increment 5vL{h). We 
have selected two scales, hi = A/1.26 and h2 — A/5.01 (A 
is the Taylor micro-scale, therefore hi is a typical scale in 
the inertial subrange and /12 is in the dissipation subrange) . 
The selection of the subensembles is performed using con- 
ditions of the type a < 5vL{h\) < b. The pdfs of 5vL{h2) 
for the unconditional ensemble and for the subensembles 
are then compared. Figure W gives the normalized pdf (cr 
is the standard deviation 01 5vL{h)) of the unconditional 
ensemble for h — hi and h — h2. We observe the clas- 
sical result that the tails of the pdfs are growing as the 
scale is decreased which is the signature of growing inter- 
mittency. The pdf also shows a negative skewness which 
is a direct consequence of the nonlinear dynamics of the 
Navier-Stokes equations. Figure ^ shows the pdfs of the 
unconditional ensemble and of a subensemble defined by 
the constraint — 1 < SvL{hi) < 0. The pdfs are almost 
superimposed, showing that the flow at scale h2 is unaf- 
fected by the condition imposed at scale hi . Figure ^ is 
a case for which the subensemble is much smaller due to 
a more restrictive condition, 1 < SvL{hi) < 4. However, 
the general behavior of the pdf supports the view that the 
chaotic dynamics of the Navier-Stokes equations tends to 
restore the original distribution. Note that the skewness is 
incorrectly predicted and seems to be correlated with the 
sign of 5vL{hi). Figure R presents a case with a very strong 
condition, —7 < 5vL{hi) < —2. Though the number of re- 
alizations is small, we observe that the top of the pdf is 
quite accurately reproduced. 

CONCLUSION 

These results, although preliminary in nature, offer cru- 
cial support to the hypothesis that a conditional average 
may be used to reduce the number of degrees of freedom 
required for the numerical simulation of turbulence. Work 
is continuing to make a more stringent assessment of the 
validity of such averages for turbulence and this includes 
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Figure 4: Normalized pdf of the unconditional ensemble. 
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Figure 5: The pdfs of the unconditional ensemble and a 
subensemble defined by the constraint —1 < SvL{hi) < 0. 
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Figure 6: The pdfs of the unconditional ensemble and a 
subensemble defined by the constraint 1 < SvL^hi) < 4. 




carrying out simulations at higher numerical resolution. At 

present we arc working on a 512'' simulation and hope to 
present results from this at the conference. 
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